/*
Este programa implementa o algoritmo de Gauss-Legendre, como
apresentado em http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm
Utiliza-se a biblioteca para aritmética de precisão GMP. A precisão mínima utilizada
em todas as variáveis relevantes do programa são tomadas como input, e interpretadas
como precisão em bits.
Realizam-se os cálculos até o valor convergir para a precisão especificada.
Além de imprimir o valor obtido em um arquivo, o programa imprime em tela
o tempo decorrido até a convergência.
*/
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

int main(int argc, const char* argv[])
{
   //Declaração de variáveis para abertura de arquivo e medida de tempo
   FILE *fp;
   char filename[50];
   struct timeval ti, tf;
   long int sec, usec;
   
   //Início da contagem de tempo
   gettimeofday(&ti, NULL);

   //Declaração das variáveis utilizadas nos cálculos
   mpf_t a, b, t, p, pi;
   mpf_t ap, bp, tp, pp;
   int precision, i;

   //Inicialização de variáveis com a precisão especificada
   precision = atoi(argv[1]);
   mpf_set_default_prec(precision);
   mpf_init(a);
   mpf_init(b);
   mpf_init(t);
   mpf_init(p);
   mpf_init(pi);
   mpf_init(ap);
   mpf_init(bp);
   mpf_init(tp);
   mpf_init(pp);
   
   //Valores iniciais da sequência
   mpf_set_ui(a, 1);
   mpf_sqrt_ui(b, 2);
   mpf_ui_div(b, 1, b);
   mpf_set_ui(t, 1);
   mpf_div_ui(t, t, 4);
   mpf_set_ui(p, 1);

   //Loop até p convergir
   for(i = 0; mpf_cmp(a, ap) != 0; i++)
   {
      //Salva os valores atuais da sequência
      mpf_set(ap, a);
      mpf_set(bp, b); 
      mpf_set(tp, t);
      mpf_set(pp, p);

      //Cálculo de a
      mpf_add(a, ap, bp);
      mpf_div_ui(a, a, 2);

      //Cálculo de b   
      mpf_mul(b, ap, bp);
      mpf_sqrt(b, b);

      //Cálculo de t    
      mpf_sub(t, ap, a);
      mpf_pow_ui(t, t, 2);
      mpf_mul(t, pp, t);
      mpf_sub(t, tp, t);

      //Cálculo de p
      mpf_mul_ui(p, pp, 2);
   }
   
   //Cálculo final de pi (pi ~= (a_n + b_n)^2/(4t_n))
   mpf_add(a, a, b);
   mpf_pow_ui(a, a, 2);
   mpf_mul_ui(pi, t, 4);
   mpf_div(pi, a, pi);
   

   //Término da contagem de tempo
   gettimeofday(&tf, NULL);
   sec = tf.tv_sec - ti.tv_sec;
   if(tf.tv_usec < ti.tv_usec)
      sec--;
   usec = (tf.tv_usec - ti.tv_usec + 1000000) % 1000000;
   printf("Tempo decorrido: %ld.%06ld s. Iterações: %d\n", sec, usec, i);
   
   //Salva o resultado para arquivo
   sprintf(filename, "pi/gauss-legendre-seq-%d", precision);
   fp = fopen(filename, "w");
   gmp_fprintf(fp, "%.Ff\n", pi);
   fclose(fp);

   return 0;
}
